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Abstract 

We present a fully grid-based approach for solving Hartree-Fock and all-electron Kohn-Sham 
equations based on low-rank approximation of three-dimensional electron orbitals. Due to the 
low-rank structure the total complexity of the algorithm depends linearly with respect to the one¬ 
dimensional grid size. Linear complexity allows for the usage of hne grids, e.g. 8192^ and, thus, 
cheap extrapolation procedure. 

We test the proposed approach on closed-shell atoms up to the argon, several molecules and 
clusters of hydrogen atoms. All tests show systematical convergence with the required accuracy. 


1. Introduction 

Electronic structure calculations arise in a variety of applications such as condensed matter 
physics or material and drug design. In this paper we focus on the development of efficient and 
accurate numerical techniques for solving Hartree-Fock (HE) [1] and Kohn-Sham (KS) equa¬ 
tions [2]. From computational point of view HE and KS equations are nonlinear eigenvalue 
problems with three-dimensional integro-differential operators. A standard way to solve these 
equations is to approximate the solution in a subspace of globally supported basis functions, e.g. 
gaussian functions. This is a classic topic and a lot of software packages are available. The 
choice of these basis functions is determined by the complexity of the the iterative process and 
the desired accuracy. It also introduces the basis set error, which might be difficult to control. 
Methods based on a sequence of embedded subspaces allow for rigorous error control. Among 
these approaches we should mention multiresolution approach based on wavelet function repre¬ 
sentation [3], hnite element methods on unstructured grids [4], finite difference method [5] and 
the projector-augmented wave method [6]. Standard hnite element or hnite difference method 
on uniform grids for a 3D-dimensional problem has Oitv' log ri) complexity (using Fast Fourier 
Transform for the evaluation of the integral operators). On non-uniform meshes we can not use 
FFT directly, thus the complexity can reach 0{rfi), where n is the one-dimensional grid size. 

One of the promising ways to reduce the complexity, proposed and studied in details in 
the papers by Khoromskij and Khoromskaia [7, 8, 9, 10, 11, 12, 13, 14, 15] is to use tensor 
decomposition of the vectors of coefficients. It was shown that in many cases, the tensor of 
coefficients can be well-approximated by the Tucker decomposition, and the storage complexity 
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goes down to 0{n log n), where n is one-dimensional grid-size. This allows for huge uniform 
grids to be used. However, to compute the solution, they still employ an intermediate basis set 
(which can be arbitrary), and successively compute the Fock matrix and update the coefficients 
with respect to this basis. The main difficulty in solving the KS/HF equation directly in the tensor 
format is that standard iterative methods should be modified in order to keep the solution in the 
tensor format, and it requires the development of a new iterative method. 

Our main contribution in this paper is an efficient 0{n log n) black-box solver which is based 
on low-rank tensor decompositions, and it does not require any additional basis sets. Thus, it 
can be used to control the basis set error in HF/KS computations. We have developed an itera¬ 
tive process that solves the original problem while all intermediate three-dimensional arrays of 
coefficients are stored in the Tucker format. Our method has several important components: we 
use the block integral iteration, fast low-rank convolution and derivative-free update of the Fock 
matrix. Finally, the extrapolation over several grid sizes is done, giving up to 0{h!^) accuracy, 
where h is the grid step. 

We organize our paper as follows. In section 2 we formulate the closed-shell Hartree-Fock 
and Kohn-Sham equations. In section 3 we discuss tensor decomposition approach and the 
Tucker format. Section 4 contains formulation of the Block-Green iteration and the derivative- 
free Fock matrix update formula. Sections 5 and 6 are concerned with obtaining the discrete 
formulation and all operations involved in the tensor formats. In section 7 we discuss the com¬ 
plexity of the proposed algorithm. Finally, in section 8 we present numerical experiments on 
closed-shell atoms up to argon, several molecules and cluster of hydrogen atoms, which illus¬ 
trate sublinear rank dependence. 

2. Hartree-Fock and Kohn-Sham equations 

Consider closed-shell electron system with Ne electrons and N — Ndl orbitals, N„uc nuclei 
with charges Zq. located at Rq., a = Then the HF/KS equations can be written as [16] 

Hmcf>i = i = UV (1) 

where 0, : ^ denotes unknown orbitals that additionally satisfy the orthogonality con¬ 

straints 

r </>i{r)(f>j(r)dr ^ 6ij. 

Jr3 

and Ai denotes orbital energies. In the Kohn-Sham model 

vm = V(p) = + Vcoulip) + VAP), (2) 


where 


N 

p(r) = 2^|0,■(r)|^ 
1=1 


is the electron density. The external potential Vext describes interaction between electrons and 
nuclei: 


Vext(r) = - ^ 
0'=1 


Zg 

|r-R.|' 


2 


(3) 



The potential VcouM is given as 


^coulix^ — 


f 

Jr3 


P(r') 

|r-r'| 


dr', 


and Vxc depends only on the density and is responsible for exchange and correlation parts of the 
operator. For the sake of simplicity we consider local density approximation (LDA) with Perdew 
and Zunger type functional [17]. We note that the concept can be extended to local spin density 
approximation (LSDA) or more accurate functionals. 

In the HF equation the potential y(<h) has the form 

vm = y„, + Vcoulip) - Km, 


where 


(pi = 


^ I *'’*',1 with T(r,r') = V^(r)0/(r'). 

Jr 3 |r - r I 

As K depends not only on density p(r), but on all the orbitals explicitly, the HF equation is 
computationally more expensive in comparison with the DFT analogues. The Fock matrix F = 
F(<1)) is dehned as 


F afi - 


- I (paHim)(ppdr, a,P-\,N. 

Js? 


It is diagonal when d) is the exact solution of (1). Total energy can be calculated as 

^ 2JJk 6 |r-r'| JJr 6 |r-r'| 

for the Hatree-Fock model and 

. . V , 1 rr p(r)p(r') _ 

^“2^/1/ oil I /I dvdv + Exci,f^ “I" Efjfi, 

2JJ^6 |r-r'| 


for the Kohn-Sham model, where 


N N 


F = 

K-inn — 




z,z. 


. , ^ IRz-R/l’ 

1=1 ]>i j 


describes the repulsion between the nuclei. 

3. Tensor decomposition approach 

Orbitals (pa, o' = 1,.. .,N depend on three variables. Discretization on a rectangular nxnxn 
mesh ' gives a three-dimensional tensor of coefficients. The storage complexity grows cubically 
in n. To get linear complexity 0(n) we will use Tucker approximation [18, 19, 20] of the all the 
tensors arising in the computations. 


'the mode sizes can be in general different, but for simplicity we always assume that they are equal to n. 
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Tensor A = {^ijkYljk^\ said to be in the Tucker format if it is represented as 


r\ r2 n 

--‘ = ZZZ Sapy^ia^ jfi^ky^ (4) 

0'=1 p=\ r=l 

where the minimal possible numbers ri, r 2 , rs required to represent A in the form (4) are called 
Tucker ranks, matrices U - V - W — {'^/tr)Sy=i called factors, 

G - {gijkYi is called core of the decomposition. 

Decomposition (4) contains r' + 3nr parameters, compared with rY elements of the whole 
tensor. In practice, instead of the decomposition we use the approximation with some accuracy e. 
Thus, each orbital on the grid is approximated in the Tucker format 

nip) nip) nip) 

a=l p=l y=l 

Our main assumption is that the Tucker ranks are small and grow only logarithmically both in 
n and e Z This has been verified before in [9]. Once the approximation is computed, we can 
check the accuracy of the computations by computing relevant physical quantities and comparing 
it with the results, obtained by other methods. 


4. Iterative method 


The standard iterative technique to solve HF/KS equations is the self-consistent field iteration 
(SCF) [16]: 


Hik) ^ik^i) ^ 


i = 1,A^. 


This is not easy to implement in the Tucker format, which is a non-linear parametrization of 
the 3D-tensor, thus the original convex minimization problem is replaced by a non-convex one. 
Flence, SCF iteration is not much simpler that the original problem. Therefore, we use a more 
convenient for our purposes Block Green iteration for Lippman-Schwinger form of (1) as it 
can be easily done within tensor arithmetics. In next subsection we describe the Block Green 
iteration technique with rotation of orbitals and present formulas to calculate the Fock matrix 
without direct Laplacian operation. 


4.1. The Block Green iteration 
Let us rewrite (1) as follows 


The action of (-A - 2/1,)“* can be written as a convolution with the Yukawa potential 

n -V^l|r-r'|| 

(-A-2/l,r'y,^(r)= —- VHr') dr' . (5) 

Jr 3 47r||r-r'|| 

We found that direct application of (5) leads to numerical difficulties, which will be discussed 
later. In this case much more efficient is to solve the screened Poisson equation directly using 
finite difference method. 
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For simplicity we start our description of the iteration for the system with one orbital (p\ = <p. 
In this case the k-th step of the Green iteration is 

0 = 2 (-A - ^ = 0/II0II- 

The energy is calculated as 

^ (6) 

The Hamiltonian = H contains the Laplacian A. The approximation error in the Tucker 
format is controlled only in the discrete L 2 norm, thus discrete differentiation will amplify the 
error. To avoid it (6) can be rewritten as 


(yd:+i)0_ 

= T® + - - - / (7) 

(0, 0) 

The generalization to the case of molecules with more than one orbital is as follows. We start by 
modifying each orbital separatedly 

= 2(-A - 2dfV‘y<*Vf^ (8) 

and orthogonalize <1) = (^ 1 ,... ,(Pn) by calculating the Cholesky decomposition of the Gram 
matrix 

r $^$£/r = LL^, 5 = (9) 

Jr3 

Then we calculate the N xN Fock matrix 


F 

= r (p^r, 

(10) 


Jk^ 


diagonalize F 

A(i+i) ^ s ’^FS, 

(11) 

and hnally rotate the orbitals by S 

,py+i) _ 

(12) 

The new values of orbital energies are 

the diagonal parts of A^*^^b 

Steps of the iterative process 


are summarized in Algorithm 1 


Algorithm 1 Block-Green iteration 
1: Calculate O = (0i,..., ^n): 0,- = 2 (-A - 
2: Orthogonalize <1): <t) = where L: LL^ = <t>^<t>dr. 

3: Calculate the Fock matrix f - Jg 3 (p dr via derivative-free formula in Statement 1 . 

4: Find new orbital energies by diagonalizing F: = S^^FS 

5: Find new orbitals: - (p^ 
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4.2. Derivative-free Fock matrix computation 

To compute the Fock matrix (10) one need to find Laplacian of orbital functions. Since all 
computations in tensor formats are done approximately, differentiation will lead to the accuracy 
loss. Therefore, we present the following derivative-free formula for the Fock matrix calculation 

Statement 1. In Algorithm 1 the Fock matrix (10) can be written in the form without derivatives 

F = f H- dr, (13) 

Jr3 ' ' 

where y® = y(<hW). 

Proof. From (10) we have 


Matrix 




Jr3 




can be written elementwise in the following way 


(14) 


Jr3 Jr3 V ^ 

Jr3 Jr3 \ ^ 

fa (v® -I- d®) (pp dx-i- A 


jf/jdx- 
^fpdx- 
® j fp dx 


Taking into account that 

^4 = (-^A-d®j y®^® 

we have 



Finally, substituting the obtained expressions into (14) we get (13). 


4.3. Density mixing 

In order to accelerate convergence we used the density mixing scheme, which is typical used 
for DFT calculations [16]. Here we provide a brief description of the scheme for the sake of com¬ 
pleteness. Our scheme is a fix point iteration which can be formally written as = G(p^^“*^). 
In the density mixing scheme, the next density p(^+0 represented as a linear combination of 
the results of previous iterations 


>=i 


1=1 
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where the coefficients a = (ai,..a^) are the solution of the minimization problem 


a = arg _min || ^ aj - G || 

7=0 

under the constraints 

7=0 


5. Discretization 


It is known that the orbitals decay exponentially at infinity [21]: 

0/(r) = , ||r|| ^ +oo. (15) 

As a result, multiplications and linear combinations of the orbitals also decay exponentially as 
||r|| ^ +O 0 . Hence, we replace the whole with a finite box Q = [-L, L]^, where L depends 
on the chosen accuracy and estimation of the smallest orbital energy We will numerically 
investigate how the choice of L in the experiments section. 

In Q. we introduce a uniform tensor-product grid lJ' = Wj' x x with h - ILilrii, where 
(Jl - {-Li + khi : k - 0 ,..., n], i = 1,2,3. We use uniform grids to obtain structured operators 
and significantly decrease computation complexity. Low-order approximation accuracy over 
uniform grids will be improved later by using extrapolation. For simplicity we use grids with 
hi - /12 = hi and Li - L 2 — L 3 . Algorithm (8)-(12) requires the computation of VO at each 
step, and the 3D-convolution have to be computed 

r 1 /^*' !ii (16) 

Jr3 l|r-r'|| 

Remark 5.1. One can use the fact that problem of finding w(r) is equivalent to find the solution 
of —Aw = 47Tf and it seems better to solve less expensive Poisson equation instead of direct 
computation of the convolution. Nevertheless, w(r) = G(l/||r||) when ||r|| —> 00 . Boundary 
conditions are unknown and, hence, to get the accuracy e one have to choose domain size L — 
0(l/e). 


On the other hand we can compute convolution with the Yukawa kernel (operator (-A - 
2/1) ' in ( 8 )) by directly inverting the shifted Laplacian. This can be done efficiently due to 
tfie exponential decay of orbitals (15), so it is safe to enforce Dirichlet boundary conditions. 
Moreover, the direct convolution with the Yukawa kernel discretized on a uniform grid leads to 
the non-symmetry of the Fock matrix and this breaks down the convergence. The point here is 
that the numerical analogue of the derivative-free formula (13) does not hold unless the action of 
(-A - 2/1) * is computed by solving screened Poisson equation. 

Thus, to discretize (-A - 2/1)“* we use standard 7-points finite difference scfieme. To dis¬ 
cretize convolutions (16) we use Galerkin scfieme and piecewise-constant basis functions fi with 
support on Qj, where le I s (0,..., n - 1 and Qi are h^ cubes centered in n = (x;, yj, Zk) e w*’. 


So, 


w(r,) = Wi^Yjfi ‘li-P 
je-T 
7 


(17) 



where 

n 0i(r)0j(r') , r 

—- —drdr, /i = I fir)(pi{r)dr. (18) 

3 l|r - r'll Jr3 

The total approximation error of the discretization scheme is O(h^). To get high-order discretiza¬ 
tion schemes translation-invariant basis functions of higher order can be used 0i(r) = (p{r - ri), 
where 0(r) is a suitable piecewise-polynomial function. 

6. Operations in the Tucker format 

In order to implement the iterative algorithm in the Tucker format, we need several basic 
operations. Linear operations (addition, multiplication by number) can be done straightforwardly 
[22]. Calculation of scalar products and norms can also be efficiently implemented within the 
tensor formats. After such operations the ranks typically increase, but can be reduced thanks to 
the SVD-based rounding procedure. 

6.1. Elementwise operations and cross approximation 

Despite linear complexity in the mode size, some operations such as elementwise products 
and convolutions may have strong rank dependence: O(r^) floating point operations. To decrease 
this complexity we use the so-called cross approximation method [23, 24, 25]. This method finds 
the decomposition using only few of elements of the approximated array. Every next iteration 
of the cross approximation procedure chooses adaptively new elements to be computed untill 
the stopping criterea holds. In the implementation we use Schur-Cross3D algorithm proposed in 
[26] that has 0{r^ + nr) complexity. 

The cross approximation algorithm is also the key technique to approximate nonlinear ele¬ 
mentwise functions of density arising in the calculation of exchange-correlation potential Vxc- 
We note, that Vxdp) itself does not have low-rank structure, however, Vxdp) 4>i is of low rank. 

6.2. Convolution 

Computation of the convolution is one of the most expensive parts of (1) with strong de¬ 
pendence on ranks of input tensors. We calculate it by the Cross-Conv algorithm proposed in 
[26]. This algorithm has 0{r* -i- nr^) complexity, which is the fastest known for practically in¬ 
teresting mode sizes n up to « ~ 2*^. However, for fine meshes it might be a good idea to use 
quantize tensor train (QTT) convolution algorithm [27] can be used. This algorithm computes 
the convolution of two tensors given in the Tucker format. The orbitals are already there, and 
it is well-known that the kernel function l/r can be approximated in the Tucker format with 
r = (9(logn log”* e). This trick with decomposition into sum of exponents was used to calculate 
external potential (3). 

6.3. Fast Laplacian inversion 

The Laplacian equation is solved using the algorithm proposed in [28]. Given the right-hand 
side in the Tucker format, the solution after the sine transform (which is separable) reduces to 
the computation of 

t^ijk 

rji + T]j + r]k-p" 

where 77 ,, / = 1,..., n are the eigenvalues of the one-dimensional Laplace operator. Such compu¬ 
tation can be straightforwardly done using the cross-approximation technique. 
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7. Complexity estimates 

For simplicity we provide an upper bound of complexity for r be the largest rank of all 
orbitals. We also assume that the calculations are done onnxnxn grid. Precomputations before 
the iterative process include calculation of the external potential (3) and the Galerkin tensor of 
the convolution kernel (18). The computation of the external potential requires N„uc ■ 0{P + 3nr). 
Precomputation of the Galerkin tensor requires Oir' + 3nr) operations. It does not depend on 
number of nuclei but has a bigger constant due to one-dimensional integration. 

Let us now estimate the complexity of one iteration. We denote by Kcross - 0{r^ -i- nr^) the 
complexity of the cross approximation calculated by Schur-Cross3D and by Kcom -0{r^ + nP' + 
rn logn) the complexity of one convolution operation [26]. We note that the constant hidden in 
0( ) is of the order of unity. Calculation of in (8) is done as follows. First of all we 

calculate Vcoui which is convolution operation that has Kcom complexity. The Coulomb potential 
is calculated only once on each iteration. In case of KS equations we run cross approximation 
procedure to calculate since the 14^ does not have low-rank structure. For the HF equa¬ 

tions we need additionally to calculate the exchange potential, which requires convolution 
operations and elementwise products. Thus, the step has NKcross complexity for KS and 
N^Kcross for the HF equations. The application of (-A - /r) ' is also of NKcom complexity but 
with approximately twice smaller constant in Kconv as one of the input tensors has fixed rank 3. 

Step (9) requires scalar product evaluations, so it has Kcross complexity as elementwise 
products are done by the cross approximation. Complexity N^Kcross is much larger than the 
Cholesky factorization as we assume that N <s: P. The Fock matrix computation (13) consists 
of which is of same complexity as yi^'O® and of scalar products. Thus, we get the overall 
complexity of one iteration to be ■ 0(P + nP + rn log ri). 

8. Numerical experiments 

The prototype is implemented in Python. The implementation of HF and KS solvers can be 
found at https; //github. com/rakhuba/tensorchem. The toolbox with arithmetics in the 
Tucker format can be found at https: //github. com/rakhuba/tucker3d. For the basic lin¬ 
ear algebra tasks the MKL library is used. Python and MKL are from the Anaconda Python Dis¬ 
tribution https : //store . continuum, io/cshop/anaconda/. Python version is 2.7.9. MKL 
version is 11.1-1. Tests were performed on 4 Intel Core i7 2.6 GHz processor with 8GB of RAM. 
However, only 2 threads were used (this is default number of threads for MKL). We would like 
to emphasize that implementation of the whole algorithm is in Python and time performance can 
be improved by implementing the most time-consuming parts of it in C or Fortran languages. 

8.1. The box size 

Here we illustrate how the results depend on the simulation box size. The orbitals decay 
exponentatially as exp(- V-2/lHOMolkll), ll-^ll ^ Therefore, one can expect that the error 
introduced by the finite size of the box has exponential decay with the box size. On Figure 8.1 
illustrates this guess. On this figure the dependence of relative error of Eh{d) with respect to the 
Eh{°°) as a function of a box size a is presented. All calculations were done for the fixed grid 
step h -Q.\ A and accuracy e - 10“®. Eh{o°) was estimated at a = 50 A. 
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Influence of a box size 



Figure 1: Dependence of refative error of £■;,(«) on a box size a for fixed ft = 0.1 A, e = 10 To calculate Ei,{oo) the 
box of size a = 50 A was used. 


8.2. Extrapolation 

As the proposed approach is fully grid-based, a good idea is to start the method on a coarse 
grid with, say A = 128 and then use it as an initial guess on finer grids (grids with 2^ grid points 
are used). Due to the linear comlexity of the proposed algorithm the total time for extrapolation 
is not larger than twice the time required for the computation on the hnest grid; 


fextrapolated - flV + tNI2 + • ■ ■ + fAf/2' - CN + CN/2 -H ■ ■ ■ -H CN< 2CN — 2tj^, (19) 


where constant C does not depend on N. This makes the extrapolation a very useful and compu¬ 
tationally efficient part of the whole algorithm. 

Numerical results showed that on coarse grids the order of the approximation is not exactly 
of the second order. Hence, we used the so-called Aitken’s delta squared process [29], which 
is equivalent to the Richardson extrapolation for the exact second order. The Aitken’s process 
accelerates the convergence of the sequence E„ such that the new sequence £' 


E'„ - E„+2 ' 


{En+2 ■ 


■En -^2 




En+2 - '2.E «+l “I" 


( 20 ) 


converges faster than £„ when n goes to inhnity. The same trick can be applied to £'. The results 
for extrapolation are shown on Figure 8.2. 


8.3. Elartree-Fock calculations on atoms 

First we present the HF calculations of closed-shell atoms He, Be, Ne, Ar. Table 8.3 repre¬ 
sents total and highest occupied molecular orbital (HOMO) energies. These results are compared 
with highly accurate results for atoms [30]. Our grid-based results are obtained with the relative 
accuracy e - 10 2. We use the relative error instead of absolute in order to obtain accurate 
HOMO values and at the same time to be faster in the calculation of lower orbitals. However, 
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Extrapolation for Ar 
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Figure 2: Extrapolation of the full energy as a function of grid size for Argon. Error is calculated with respect to highly 
accurate results from [30]. 


6 

10-" 

10-" 

10-' 

10-"' 

10-" 

Total energy 

-2.8615 

-2.861678 

-2.8616801 

-2.861680000 

-2.86167999593 


Table 1: Helium total energy for different values of relative error e. The energy was extrapolated on a sequence of grids: 
from 128^ to 8196^. Bold denotes correct numbers. The total energy in [30] is -2.861 679 996. 


absolute accuracy can also be used. Extrapolation is done on a range of grid sizes: from 128^ to 

81923 . 

Results in Table 8.3 show systematic convergence of the total and HOMO energies. Note, 
that for all atoms under consideration extrapolation with grid sizes up to = 8192^ grid points 
is enough to get accuracy e - lO”^. Without extrapolation, much larger grids are needed. 

8.4. LDA calculations on molecules 

We perform the LDA calculations for several molecules. Molecular geometries were taken 
from the NIST database [32]. In Table 8.4 the resulting total and HOMO energies are compared 
with the LDA computations with the aug-cc-pVXZ (X=Q, 5) basis sets. The latter we done using 
NWCHEM program package [31]. In all experiments the relative error for all orbital functions 
was set to be e = 10 3. The box size was chosen adaptively to the HOMO energy. 

We note that the extrapolated energy is smaller than for the basis aug-cc-pV5Z approximately 
10 hartree both for the total energy and for the HOMO energy, thus our method has less basis 
set error. Eigure 8.4 illustrates rapid convergence of all hve occupied orbitals with respect to the 
number of iterations for e = 10 3. As anticipated the convergence stopped below e. We note that 
in all experiments orbitals have converged to the e = 10 ^ within 30 iterations. 
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Atom 

Method 

Total energy 

HOMO energy 

He 

N = 8192 

-2.861 670 

-0.917 950 


Extrapolated 

-2.861 680 

-0.917 955 


[30] 

-2.861 680 

-0.917 956 


aug-cc-pVQZ 

-2.861 539 

-0.917 915 


aug-cc-pV5Z 

-2.861 635 

-0.917 935 

Be 

N = 8192 

-14.572 256 

-0.309 263 


Extrapolated 

-14.573 023 

-0.309 269 


[30] 

-14.573 023 

-0.309 270 


aug-cc-pVQZ 

-14.572 976 

-0.309 269 


aug-cc-pV5Z 

-14.573 011 

-0.309 270 

Ne 

N = 8192 

-128.518 74 

-0.850 523 


Extrapolated 

-128.547 08 

-0.850 410 


[30] 

-128.547 09 

-0.850 410 


aug-cc-pVQZ 

-128.544 69 

-0.850 210 


aug-cc-pV5Z 

-128.546 87 

-0.850 391 

Ar 

N = 8192 

-526.740 5 

-0.591 024 


Extrapolated 

-526.817 4 

-0.590 017 


[30] 

-526.817 5 

-0.591 017 


aug-cc-pVQZ 

-526.816 9 

-0.591 013 


aug-cc-pV5Z 

-526.817 3 

-0.591 011 


Table 2: Total and HOMO energies for dilferent atoms, 6 = 10 Extrapolated calculations were done on a sequence of 
grids: from 128^ to 8196^. aug-cc-pVQZ and aug-cc-pV5Z calculations were done by NWCHEM [31]. 


Molecule 

Method 

Total energy 

HOMO energy 

H 2 

N = 8192 
Extrapolated 
aug-cc-pVQZ 
aug-cc-pV5Z 

-1.137 392 3 
-1.137 392 8 
-1.137 249 9 
-1.137 374 8 

-0.378 667 
-0.378 668 
-0.378 649 
-0.378 665 

CH 4 

N = 8192 
Extrapolated 
aug-cc-pVQZ 
aug-cc-pV5Z 

-40.116 452 
-40.119 829 
-40.118 644 
-40.119 299 

-0.348 989 
-0.348 984 
-0.348 964 
-0.348 982 

C 2 H 6 

N = 8192 
Extrapolated 
aug-cc-pVQZ 
aug-cc-pV5Z 

-79.069 964 
-79.075 142 
-79.070 784 
-79.072 763 

-0.299 763 
-0.299 761 
-0.299 724 
-0.299 762 


Table 3: Total and HOMO energies for different molecules, relative error 6 = 10 Extrapolated calculations were done 
on a sequence of grids: from 128^ to 8196^. aug-cc-pVQZ and aug-cc-pV5Z calculations were done by NWCHEM [31]. 
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Figure 3: Relative errors of each orbital energy as a function of number of iterations (e = 10 ^). 


8.5. LDA calculations on cluster of atoms 

Here we present calculations on a test system of hydrogen atoms, which was firstly calculated 
in [13] as an artificial example. The system represents finite cluster of hydrogen atoms with 
unit cell having primitive cubic structure. The distance between atoms is set to be tf = 2 A. 
As in section with LDA calculations on molecules in all experiments we observed systematical 
convergence of the total energy with e precision. Table 8.5 illustrates that ranks of orbitals grow 
sublinearly with the system size. Thus, we expect the proposed algorithm to be particularly 
efficient for systems with regular location of atoms. The time of one iteration on A = 1024 and 
e - 10'^ is 9 sec for H^xixi and 68 sec for Hgxixi, which as anticipated scales approximately 
quadratically with the size of the system: (8/3)^ » 7.1 and 68/9 » 7.5. 

9. Conclusions and future work 

Tensor approach to the solution of 3D problems in electronic computations allows us to reach 
desirable accuracy at moderate computational cost. This can be used, for example, in verification 
of other methods and in construction of accurate basis sets for other methods [33]. The main 
problem with the presented approach is still in its scaling for large N: for example, one has to 
compute the approximation of all the products fifj using the cross method. This can be solved 
by the resolution of identity methods [34], but they are typically used for global basis sets. What 
is interesting, is that the density itself shows a very good low-rank structure; thus, it is very 
interesting to apply tensor decomposition methods for the orbital-free DFT functionals [35, 36]. 
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Cluster 

min ranks orbital 

max ranks orbital 

Hix2x2 

16 X 16x 14 

17 X 17x 15 

H3x2x2 

28x28x21 

35 X 35 X 20 

H8x2x2 

25 X 25 X 22 

36 X 36 X 26 

Hix2x1 

13x 12x 12 

13 X 12x 12 

H2x2x1 

16 X 16x 14 

17 X 17x 15 

H5x2x1 

16x 18x 11 

16 X 20 X 13 

H9x2x1 

14x 18x 11 

17 X 22 X 13 

Hi6x2x1 

16x 19x 11 

21 x27x 14 

Hix4x1 

13 X 12x 13 

13 X 13 X 13 

H3x4x1 

19 X 20 X 13 

23 X 23 X 15 

H8x4x1 

23 X 20 X 12 

32 X 27 X 15 


Table 4: Ranks of orbitals with maximum and minimum ranks on different clusters of hydrogen atoms. 
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